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We present the results of a molecular dynamics simulation study of thermodynamic and structural 
properties upon supercooling of a low concentration sodium chloride solution in TIP4P water and 
the comparison with the corresponding bulk quantities. We study the isotherms and the isochores 
for both the aqueous solution and bulk water. The comparison of the phase diagrams shows that 
thermodynamic properties of the solution are not merely shifted with respect to the bulk. Moreover, 
from the analysis of the thermodynamic curves, both the spinodal line and the temperatures of 
maximum density curve can be calculated. The spinodal line appears not to be influenced by the 
presence of ions at the chosen concentration, while the temperatures of maximum density curve 
displays both a mild shift in temperature and a shape modification with respect to bulk. Signatures 
of the presence of a liquid-liquid critical point are found in the aqueous solution. By analysing 
the water-ion radial distribution functions of the aqueous solution we observe that upon changing 
density, structural modifications appear close to the spinodal. For low temperatures additional 
modifications appear also for densities close to that corresponding to a low density configurational 
energy minimum. 

PACS numbers: 61.20. Ja,64.60.My,61.20.-p 



I. INTRODUCTION 

The anomalous behaviour of water has attracted the 
attention of scientists since long time. 

Many fascinating features distinguish water from most 
other liquids. Among them we point out the presence 
of a line of temperatures of maximum density (TMD) in 
the P — T plane and the power law divergences upon su- 
percooling of the coefficient of thermal expansion, of the 
specific heat and of the isothermal compressibility^! 3 -! 4 - 

Molecular dynamics simulations of water modeled with 
ST2, SPC/E, TIP4P and TIP5P potentials have shown 
the existence of a second critical point^&^&^ i 10 ' 11 buried 
inside the supercooled region with a coexistence between 
a low density liquid (LDL) and a high density liquid 
(HDL). The deep supercooled region is not directly acces- 
sible to experiments, nonetheless signatures of the second 
critical point have been observed* 1 ^ In this picture wa- 
ter anomalies are a direct consequence of the existence 
of the second critical point in the region of low tempera- 
tures and high pressures, the water spinodal line is non- 
reentrant and correspondingly the TMD line bends to 
avoid crossing the spinodal line. 

Being water a rather unique solvent, the properties of 
aqueous solutions are also a topic of substantive interest. 
In particular thermodynamic and structural properties 
of aqueous solutions of electrolytes are of importance be- 
cause of the ubiquitous presence of these solutions in bi- 
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ological and geophysical systems. Besides the properties 
of aqueous solutions upon supercooling are significant for 
cryopreservation of biological tissues. 13 ! 14 

As far as water itself is concerned it is well known that 
the presence of solutes affects important properties like 
melting point, boiling point and viscosity. It is therefore 
important to understand to what extent thermodynamic 
properties of water are affected by the solutes upon su- 
percooling. In fact while many studies in literature con- 
centrate on ionic solutions at ambient temperature, very 
little has been done in the supercooled region. Differ- 
ential scanning calorimeter measurements have shown 1 ^ 
that for low NaCl content, up to about 1 mol/kg, at am- 
bient pressure, the TMD and the anomalous behaviour 
of specific heat upon supercooling are still present, grad- 
ually tending to disappear upon increasing the salt con- 
tent. 

A large number of experiments and simulations focus 
on the hydration structurei^^^^^^^^ 3 .^^^^ at 
ambient temperature of ionic solutions, that is to say on 
the arrangement of water molecules in proximity to the 
ions. At low concentrations the presence of the ions per- 
turbs the structure of water only locally, while it remains 
unaltered far away from the ions^ 

Traditionally, ions have been divided in structure- 
making (or cosmotrope) and structure-breaking (or 
chaotrope) depending on the effect they have on water- 
water hydrogen bonds. The first appearance of these 
concept dates back to the Hofmeister's series^ More 
than thirty effects connected to the Hofmeister's series^ 
can be found, nevertheless its greatest significance resides 
in the different degree of hydration of ionic compounds 
in aqueous solutions. Although some authors 3 ^ ques- 
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tioned the significance of the structure-making/breaking 
concept, mainly because of the difficulty in quantifying 
correlated properties, the majority of studies apply the 
usual definitional Structure-making ions break the hy- 
drogen bonds of surrounding water molecules, thus al- 
lowing a rearrangement of water molecules that results in 
a ordered hydration structure. Structure-breaking ions, 
bend without breaking the hydrogen bonds of nearby wa- 
ter molecules, thus disordering the tetrahedrally coordi- 
nated network. So far it is not clear whether and how 
structure-making/breaking properties of ions are con- 
nected to modifications of thermodynamic properties of 
the solvent upon supercooling. 

We present here a molecular dynamics (MD) simula- 
tion study of an aqueous sodium chloride solution, in 
the following denoted as NaCl(aq), upon supercooling. 
The salt content in the aqueous solution amounts to 
c = 0.67 mol/ kg. Its thermodynamic and structural 
properties are compared to a MD simulation of bulk wa- 
ter. The isotherms and the isochores are studied for both 
systems to investigate whether and to what extent the 
thermodynamic properties of water upon supercooling 
are modified in presence of a low concentration of ions. 
The analysis of the phase diagrams leads to the determi- 
nation the TMD line and the spinodal line. 

The structural properties are also studied, calculating 
the partial radial distribution functions, RDFs, of wa- 
ter molecules around the ions. The RDFs features are 
analysed also taking into account the structure-making 
properties of the Na + and the structure-breaking prop- 
erties of the Cl~ ions. 

The paper is organized as follows: in Sec.[ll]we report 
the simulation details. The thermodynamic results are 
presented in Sec. lIIII and discussed therein. In Sec.lIVIwe 
present the hydration radial distribution functions and 
we discuss their trends as function of density. Conclu- 
sions are drawn in Sec. W\ 



II. MODEL AND COMPUTER SIMULATION 
DETAILS 

Molecular dynamics simulations are performed on two 
systems, bulk water and water with Na + and Cl~ ions. 
The systems are simulated with the DL_POLY package^! 
using the TIP4P site model for water molecules^ 

The particles interact via short range L-J potentials 
and long-range coulombic potentials. Thus the interac- 
tion potential can be written as: 



parameter e,j and are obtained using the Lorentz- 



Bcrthclot mixing rules e»j = 



Ej and Oij 



V, l{ r) = if 



4e i: 
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where ion-water and ion-ion parameters are taken from 
Koneshan and Rasaiah calculation for L-J potential, 16 
based upon a re-parametrization of Pettitt and Rossky 
parameters^ for the Huggins-Mayer potential. The L-J 



The L-J interaction parameters and the charges are 
shown in Table |T] along with masses and charges of the 
sites. 



TABLE I: Masses, electric charges and L-J interaction param- 
eters for TIP4P water sites O, H and X and ions. 








H 


X 


Na+ 


cr 


m (a.m.u) 


16.000 


1.008 




35.453 


22.990 


q(e) 




+0.52 


-1.04 


+1.00 


-1.00 


eo-j (kJ /mol) 


0.65000 






0.56014 


1.50575 




3.154 






2.720 


3.550 


tH-j(kJ /mol) 








0.56014 


1.50575 


a H -j{k) 








1.310 


2.140 


e Na+-j( fcJ / mo 


0.56014 


0.56014 




0.11913 


0.35260 


°Na+-j(A) 


2.720 


1.310 




2.443 


2.796 


e cl - _-{kJ /mol) 


1.50575 


1.50575 




0.35260 


0.97906 




3.550 


2.140 




2.796 


3.487 



The cut-off radius is set to 9.0 A and the long range 
coulombic interactions are treated by the Ewald summa- 
tion method, with the convergence parameter oM set to 
6.4/L, where L is the edge length of the cubic simulation 
box in which the particles reside. 

The MD simulations are carried out in the NVT ensem- 
ble. The temperature is fixed by the use of the Berendsen 
thermostat.— The integration timestep employed is 1 fs. 

In the case of bulk water the system is simulated 
with 256 water molecules while in the NaCl(aq) system 
we have 250 water molecules plus 3 sodium ions and 
3 chloride ions. The corresponding concentration ex- 
pressed as moles of solute per mass of solvent (molality) 
is c = 0.67 mol /kg. 

We span densities in the range 0.80 g/cm 3 < p < 
1.1 g/cm 3 for the NaCl(aq) and 0.80 g/cm 3 < p < 
1.05 g/cm 3 for bulk water. The analysed densities and 
the corresponding boxlenghts are reported in Table |TT] 
for both systems. In the starting configurations particles 
lie on a face centered cubic lattice with random orienta- 
tion of water molecules. The crystal is melted at 500 K 
and temperature is then progressively decreased during 
the equilibration runs. The longest equilibrations have 
been run for 10 ns. For both systems we examine tem- 
peratures in the range 210 K < T < 500 K. Some of 
the thermodynamic points for the bulk system are taken 
from a previous work of our group i^ 6 . 



III. THERMODYNAMIC RESULTS 

We investigate how the thermodynamic behaviour of 
water upon supercooling changes in the presence of a low 
concentration of Na + and Cl~ ions with respect to the 
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TABLE II: Densities and boxlenghts for NaCl(aq) and bulk 
water. 



n ( o /cm?} 


NaClfaa) L h ^(k) 




1.1 


19.189 




1.05 


19.489 


19.388 


1.025 


19.641 




l.U 




i n 71 n 

iy. 1 1U 


0.98 


19.943 


19.839 


0.95 


20.151 


20.046 


0.90 


20.517 


20.410 


0.87 


20.750 


20.650 


0.85 


20.911 


20.811 


0.83 




20.977 


0.80 


21.338 


21.236 



bulk by looking at isotherms and isochores obtained from 
simulated state points. 

From the isotherms the limit of mechanical stability 
can be determined in the framework of mean field theo- 
ries from the divergence of the isothermal compressibility. 
The singularity points where 



dP 
~8p~ 



= 



(2) 



define the spinodal line. They can be obtained searching 
for the minima of the Pt{p) curves. 



In Fig. Q] we show the isotherms of the NaCl(aq) sys- 
tem obtained upon cooling for densities ranging from 
p = 0.80 g/cm 3 to p = 1.1 g/cm 3 and temperatures 
from T = 300 K to T = 210 K. The spinodal line 
has been drawn joining the minima of the isotherms. 
The isotherms from T = 500 K to T = 350 K do not 
show a minimum and they have been not reported. The 
T = 300 K isotherm is the only one to have the minimum 
in correspondence of the density p = 0.85 g/cm 3 , all the 
other lower temperature isotherms show a minimum cor- 
responding to the density p — 0.87 g/cm 3 . 

We report in Fig. [2] the isotherms of bulk water in 
the P — p plane. We plot the bulk isotherms for densi- 
ties ranging from p — 0.80 g/cm 3 to p = 1.05 g/cm 3 and 
temperatures ranging from T — 300 K to T = 210 K. 
Analogous to what happens for the solution the highest 
temperatures isotherms, T = 500 K, T = 400 K and 
T = 350 K do not exhibit minima (not shown). The 
spinodal line for bulk water, drawn joining the isotherms 
minima, shifts from p — 0.83 g/cm 3 , for T ranging from 
300 K to 260 K, to p = 0.85 g/cm 3 , for the lowest tem- 
peratures. 

The comparison of the region of the P— p plane close to 
the spinodal line of the two systems is reported in Fig. [3] 
The density range is here 0.80 g/cm 3 < p < 0.95 g/cm 3 . 
We observe that the NaCl(aq) isotherms minima are 
shifted to higher densities with respect to the bulk, more- 
over the curves envelope is globally shifted to lower pres- 
sures. We note that the presence of the ions also changes 
the shape of the curves. The isotherms of the NaCl(aq) 
system appear in fact closer to each other than in the 
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FIG. 1: Spinodal line and isotherms of the NaCl(aq) solution FIG. 2: 
in the P — p plane. plane. 



Spinodal line and isotherms of bulk water in the P — p 



4 



bulk phase. 

In Fig. H] we present a zoom of the region of the P — p 
phase diagram in the density range 0.90 g/cm 3 < p < 
1.05 g/cm 3 , in order to compare the isotherms of the two 
systems at high densities. In both systems the T = 220 K 
and the T = 210 K curves cross the higher temperature 
isotherms and change the sign of their curvature. These 
effects are particularly evident for the bulk T = 220 K 
and T = 210 K isotherms and the NaCl(aq) T — 220 K 
isotherm. In bulk water these changes of curvatures have 
been already observed^i and related to the existence of 
a second critical point. In fact they are signatures of the 
approach of the system to the coexistence between low 
an high density liquid phases. Therefore we hypothesize 
also for the NaCl(aq) the existence of a second critical 
point. For the change of curvature a shift toward higher 
densities can be noticed in the aqueous solution. This 
observation leads to infer that also the second critical 
point might consequently be shifted to a higher density 
with respect to the position in the bulk. 
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FIG. 4: Blow up and comparison of the NaCl(aq) and 
bulk isotherms in the P — p plane, in the density range 
0.90 g/cm 3 < p < 1.05 g/cm 3 . Symbols are as in Fig. [T] and 
Fig.H 
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FIG. 3: Blow up and comparison of the NaCl(aq) and bulk 
spinodal line and isotherms in the P — p plane, in the density 
range 0.80 g/cm 3 < p < 0.95 g/cm 3 . Symbols are as in Fig.fl] 
and Fig. [2] 

We now consider the isochores of the NaCl(aq) and the 
pure water systems reported in Fig. and Fig. [S] respec- 
tively. The behaviour of the TMD line can be extracted 
from the analysis of the isochores plane. In fact along 
the TMD line the coefficient of thermal expansion 



ap 



1 fdV 



V V dT 



goes to zero. The thermal pressure coefficient 



dP 
df 



(3) 



(4) 



iv 



ap 



(5) 



is connected to the coefficient of thermal expansion by 



where Kt is the isothermal compressibility. Therefore 
the TMD points lie on the line connecting the minima of 
the isochores. For both systems the minima have been 
determined fitting the isochores with a fourth order poly- 
nomial function. 

For the aqueous solution the isochores in Fig. \5\ are 
shown from p = X.lg/cm to p = 0.87 g/cm 3 in the 
temperature ranges from T = 500 K to T — 210 K. 
The p = 0.87 g/cm 3 isochore almost coincides with the 
spinodal line (see Fig. Q]) also reported in Fig. [5] In 
Fig. [S] the isochores of the bulk water are reported from 
p = 1.05 g/cm 3 to p = 0.85 g/cm 3 , in the same temper- 
ature range as above. Here the p = 0.85 g/cm 3 isochore 
almost coincides with the water spinodal shown in the 
same figure. For both systems the isochores with den- 
sities lower than that of the the spinodal lines are not 
reported since they lie beyond the mechanical stability 
limit of the system. 

For bulk water we observe a spinodal extending to the 
region of negative pressures and non reentrant down to 
the lowest temperature investigated. This behaviour has 
been already observed in literature for different poten- 
tials J 6>8 i 37 1 3S i 3g 1 40 1 4i 1 42 AlgQ ^ aqueous so l u tion displays 

a non reentrant spinodal line qualitatively similar to the 
bulk. 

In the case of the aqueous solution all the isochores 
above p = 0.87 g/cm 3 exhibit a clear minimum. This 
result agrees with the experimental work by Archer and 
Carter— in which they find a TMD line for the solution 
at the concentration studied here. From the minima, via 
a fitting procedure, we determine the exact location of 
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FIG. 5: Spinodal line, TMD line and isochores of the 
NaCl(aq) solution in the P — T plane. 



FIG. 6: Spinodal line, TMD line and isochores of bulk water 
in the P — T plane. 



the TMD. We observe that the TMD line, reported in 
Fig. [5j reverts its path avoiding the spinodal. 

In the case of pure water (Fig. [6]) only for the first 
four isochores a minimum is found in the range of tem- 
peratures explored. The shape of the TMD line is in 
agreement with the results that can be found in litera- 
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The comparison of the isochores planes of NaCl(aq) 
(Fig. [5]) and bulk water (Fig. [6} reveals that in the P — T 
plane the isochores of NaCl(aq) lie at lower pressures 
than the corresponding bulk curves. Besides, the concav- 
ity of the isochores in proximity to the minima is deeper 
in bulk water than in NaCl(aq). This is mirrored by the 
different pattern followed by the TMD line in the two 
systems. 

The location of the second critical point for bulk 
TIP4P water is estimated to be at negative pressure by 
Tanaka 6 , while it is found at positive pressures for ST3& 
TIP5P— and TIP4P-ewii. From the behaviour of our is- 
cohores, TMD and spinodal for both the bulk water and 
the aqueous solution it is actually not easy to establish 
whether it falls at positive or negative pressures. 

In Fig.[7]the TMD line of NaCl(aq) and bulk water are 
reported together with the spinodal lines of the two sys- 
tems. Comparing the TMD lines we notice that, apart 
for the mild temperature shift due to the colligative prop- 
erties of the solutions, the TMD lines appear different in 
shape. In fact, spanning the same range of densities, the 
NaCl(aq) TMD extends in a region of pressures wider 
than that of the bulk. Therefore, the effect of the ions 
upon supercooling is not limited to give a rigid transla- 
tion of the phase diagram. 
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FIG. 7: TMD and spinodal lines for NaCl(aq) and bulk water 
in the P — T plane. 



While the presence of the ions seems to induce impor- 
tant changes in an anomalous thermodynamic property 
of water, the TMD, the limit of mechanical stability ap- 
pears not affected by the ions. In fact from Fig. [7]we can 
see that the spinodal lines of the two systems are almost 
indistinguishable from each other. 
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In Fig. [8] we present the configurational (or potential) 
energy (per molecule) of the NaCl(aq) and bulk sys- 
tems, at T = 220 K, as a function of density. In the 
NaCl(aq) U presents a minimum in correspondence of 
p = 0.95 gj err? while in the bulk the minimum is reached 
for p = 0.86 g/cm 3 . These energy extrema can be related 
to the existence of LDL and HDL phases. It has been 
shown by Kumar et that two minima are present 
for TIP5P water confined between hydrophobic plates for 
p = 0.88 g/cm 3 and for p = 1.39 g/cm 3 . The presence of 
this minimum in our system supports the existence of a 
liquid-liquid (LL) coexistence also in the aqueous solu- 
tion. 
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FIG. 8: Comparison of the configurational energy per 
molecule, at T = 220 K, as a function of the density. 



IV. STRUCTURAL RESULTS 

We discuss in this section the radial distribution func- 
tions (RDFs) of the aqueous solution, g(r), computed 
averaging the equilibrated configurations. 

In the NaCl(aq) system at the salt content investi- 
gated we do not observe significant changes in water- 
water RDFs with respect to the bulk system. Minor dif- 
ferences occur only in the secondary peaks, for low den- 
sities and in the deep supercooled region (not shown) . It 
has been shown in a number of paper o 17 i 18 ' 19 i 20 that the 
effect of the ions on water-water structure is similar to 
the one obtained applying an external pressure and that 
it becomes more evident upon increasing the salt content. 
The salt concentration in our system, c = 0.67 mol/kg, is 
probably too low to find clear evidence of this effect. Be- 
sides both sodium and chloride ions can be found in the 
central region of the Hofmeister series so that not very 
marked effect on water-water structure are expected. 

We present in the following the analysis of the water- 
ion structure comparing RDFs along the isotherms T — 



300 K and T = 220 K for densities from p = 1.1 g/cm 3 to 
p = 0.80 g/cm 3 . We note that for the highest density at 
T = 300 K the RDFs show different trends with respect 
to all other densities investigated. These trends likely 
anticipate a different behaviour of the system at very 
high densities (not investigated in the present paper). 

In Fig. [9] and Fig. [10] we show the sodium-oxygen and 
the sodium- hydrogen RDFs. 

We observe that for both temperatures the first and 
the second coordination shells of the gNa-oir) are well 
defined showing the tendency of water to form stable 
hydration shells around the sodium ion, usually classified 
as a structure-making. For the gNa-H{r) the first shell 
is also well defined while the second shell is broader. 

The position of the first peak of the gNa-H{ r ) is shifted 
to a higher distance with respect to the position of the 
first peak of the gNa-o( r )- The difference between the 
two peaks is circa 0.6 A indicating that the two hydrogens 
of the molecule lie symmetrically around the sodium ion. 
This result is confirmed by the coordination numbers of 
the first shell which we report in Fig. [TT] together with 
the coordination numbers of the Cl-0 and Cl-H. 

We consider now more in detail the behaviour of the 
9Na-o( r ) a s function of the density. 

At T = 300 K starting from p = 1.05 g/cm 3 the peaks 
heights decrease significantly until they reach a minimum 
at the density of the spinodal line, p — 0.87 g/cm 3 . For 
densities lower than this value down to p = 0.80 g/cm 3 
the peaks heights slightly increase again. On going 
from high density down to the spinodal density along 
the isotherm the Na + ions partially lose their structure- 
making ability due to the approach to the instability re- 
gion, where the onset of large fluctuations and a tendency 
to demixing could appear. 

At T — 220 K the two peaks slightly sharpen as a 
consequence of decreasing temperature. Starting from 
p = 1.05 g/cm 3 the height of the first peak of the 
9Na-o( r ) now shows a maximum in proximity of the en- 
ergy minimum density, then it decreases showing a min- 
imum on approaching the spinodal density. We observe 
close to the second peak the appearance of secondary 
shells especially for densities close to the spinodal and 
to the energy minimum density, e.g. a possible region of 
LL coexistence. This indicates that the second shell of 
hydration is less defined. 

The gNa-nif), shown in Fig. 1101 displays an analogous 
trend as a function of density both for T = 300 K and 
T = 220 K with respect to that observed in Fig. [5] for the 
sodium-oxygen pair. In particular we observe also here 
the appearance of secondary shells. 

In Fig. [T2l and Fig. [13] we show the gci-o{ r ) and the 
gci-u(r). 

For the gci-o( r ) we see a broad and structured first 
shell and a not well defined second shell. For the 
9ci-h{i") we observe a well defined first peak and a 
broad and structured second shell. At variance with the 
9Na-o( r ) an d the gNa-H(r) the first and the second shell 
are not always well defined for these RDFs because of the 
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FIG. 9: Na-O radial distribution functions at T = 300 K (top 
panel) and T = 220 K (bottom panel). The curves are shifted 
by a unit step on the y-axis for the sake of clarity. 



structure breaking property of the CI ion. 

By comparing the positions of the first peaks of the 
gci-o( r ) and the gci-H(r) we note an approximated 
charge ordering around the negative chloride ion. The 
partially positive charged hydrogens move closer to the 
ions and the first shell of the oxygens is placed around 
the first minimum of the gci-H(r)- 

We now analyse the behaviour of the Cl~ RDFs as 
a function of densities for the two temperatures investi- 
gated. 

For the gci-o(r) at T = 300 K the first shell of oxy- 
gens is well defined only at p — 1.05 g/cm 3 . Below 
p = 1.05 g/cm a shoulder appears on the right side of 
the first peak. The peaks heights of the gci-o( r ) show a 
trend similar to that of the gNa-o(f) while their shape 
is more strongly dependent on the density. On approach- 
ing the spinodal the first peak merges with the shoulder 
resulting in a broad structure. The structure-breaking 
role of the Cl~ ions seems to be reinforced close to the 
instability limit. 

At T — 220 K the shoulder close to the first peak 
evolves in a second peak indicating a less stable configu- 
ration of the chloride hydration with respect to sodium. 
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FIG. 10: Na-H radial distribution functions at T = 300 K 
(top panel) and T = 220 K (bottom panel). The curves are 
shifted by a step of 0.5 on the y-axis for the sake of clarity. 



In fact the first hydration shell is characterized by the ex- 
istence of two possible equilibrium positions albeit close 
to each other for the chloride-oxygen pair. In this first 
shell oxygens of two kinds are present, those tied to the 
hydrogens of the first shell of the gci-H and those tied 
to hydrogens belonging to the second shell of the gci-H 
that in fact also shows double peak structure (see Fig fTB")) . 
For this temperature a structuring of the curves also ap- 
pears for densities close to the energy minimum density 
and to the spinodal density. 

For the gci-H{r), Fig. [131 the trend of the peaks in- 
tensities as a function of density are similar to those in 
Fig. H2] but the RDFs appear much more structured with 
a second and third shell evident at all densities. The 
second shell is modulated in a double peak structure for 
both temperatures. We also note that for both tempera- 
tures the first shell of hydrogens of the gci-H{r) is well 
defined at all densities, at variance with the second peak. 
It shows that hydrogens penetrate inside the shell of the 
oxygens for the attraction of the ions. 

The broadening and the modulation of the oxygen 
first peaks and the hydrogen second peaks in the the 
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FIG. 11: Coordination numbers relative to the first shell for 
the Na-O, Na-H, Cl-0 and Cl-H pairs. 



9ci-o{ r ) an d the gci-H( r ) respectively seem to be re- 
lated to a competition between the formation of the hy- 
drogen bonds and the repulsion/ attraction of the chloride 
ion with respect to oxygen/hydrogen atoms. 

We have seen that the Na + ions are able to form a suf- 
ficiently rigid shell of oxygens, see the first peaks in Fig. [SI 
breaking the hydrogen bonds when necessary while the 
Cl~ ions can attract the hydrogens in a first shell and 
repel the oxygens, but they have little effect on the hy- 
drogen bonding. Coordination numbers confirm this pic- 
ture, see Fig. 111! They show in fact no correlation in 
trend between the first shell of hydrogen and the first 
shell of oxygens around the chloride ion as a function of 
density. 

For each of the four hydration RDFs the peaks 
positions arc compatible with what found in litera- 

The ion-ion RDFs (not shown) have a low statistic due 
to the low concentration. The positions of the peaks 
indicate the presence of clusters of positive and nega- 
tive ions in agreement with what reported in several pa- 



pers. 



16,25,26,46 



both theoretical and experimental. 



V. CONCLUSIONS 

We performed molecular dynamics simulations of a 
sodium chloride aqueous solution with concentration c = 
0.67 mol/kg and of bulk water upon supercooling at dif- 
ferent densities. We used TIP4P water model both for 
the aqueous solution and for the bulk. 

We studied how the thermodynamic behaviour of bulk 
water is modified by the presence of ions and in particular 
the possible modification of the TMD and of the spinodal 
line. We found that while the presence of ions influences 
the shape and the position of the TMD, it leaves sub- 
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FIG. 12: Cl-O radial distribution functions at T = 300 K 
(top panel) and T — 220 K (bottom panel) The curves are 
shifted by a unit step on the y-axis for the sake of clarity. 



stantially unaltered the spinodal line with respect to the 
bulk. Besides the TMD curve appears broadened and 
it extends to much lower temperatures than bulk TMD. 
In a previous study of TIP4P water under hydrophobic 
confinement^ it was found that both the spinodal line 
and the TMD are shifted toward higher pressures. In the 
case of the aqueous solution studied in the present work 
the TMD line shifts also to lower temperatures and corre- 
spondingly its curvature broadens. In addition, the pres- 
ence of the ions moves the isotherms and the isochores to- 
ward lower pressures and it slightly modifies their shape, 
especially in proximity to the minima of those curves. 
Therefore, although the limit of mechanical stability of 
the system remains unaltered in the solution, the effect 
of ions is not limited to a trivial pressure shift. The 
modification of the isotherms and the isochores curves 
and especially of the TMD shows that in the supercooled 
regime the ions affect the thermodynamics of the system 
even at this low concentration. 

In our aqueous solution we observed at higher densi- 
ties a crossing of the two lowest temperatures isotherms 
with the higher temperatures ones. It has been shown 
that the change of sign of the curvature of isotherms in 
this density region is connected to the existence of a sec- 
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FIG. 13: Cl-H radial distribution functions at T = 300 K 
(top panel) and T = 220 K (bottom panel). The curves are 
shifted by a step of 0.5 on the y-axis for the sake of clarity. 



ond critical point in bulk water M&L We can consequently 
hypothesize the existence of a LL critical point also in 
the solution probably shifted to higher densities. The 
TMD shape is related to the position of the LL critical 
point hence our results indicate that the presence of ions 
stabilizes the high density liquid phase driving down in 
temperature the LL line. Interestingly two recent exper- 
imental papers by O. Mishim a 47 ' 48 on LiCl- water solu- 
tions seem to indicate a suppression of the transition to 
the low density amorphous-like (LDA-like) water. This 
appears consistent with our TMD line shift. Besides, 
according to J. Holzmann et al. the shift to low temper- 
atures of the TMD can be framed in a picture in which 
the free water, that is the water far from ions, becomes 
upon supercooling more mobile and pressurizedi 4 ^ 

The spinodal line, related to low density liquid prop- 
erties, does not show any substantial modification with 
respect to the bulk at least when a small amount of ions 
is present. 

The shape of the curve describing the configurational 



energy as a function of the density (see Fig. [8]) for the 
NaCl(aq) closely resembles what found by Kumar et al^ 
for confined water in the low density region zone explored 
in this work. Importantly we observe the presence of a 
minimum analogous to the LL low density coexistence 
minimum observed in bulk water. 

The structural analysis of the NaCl(aq), shown here 
for T = 300 K and at T = 220 K, revealed that the 
four hydration radial distribution functions gjq a -o{ r )- 
9Na-H{r), g C i-o{r) and gci-H{r) display minima in 
peaks intensities in proximity to the density correspond- 
ing to the spinodal line at both high and low tempera- 
tures and the appearance of several secondary peaks for 
densities close to the spinodal line and to the configura- 
tional energy minimum, at low temperatures. Therefore 
we found that the spinodal line and the energy mini- 
mum, low density LL coexistence, affect significantly the 
hydration structure. 

As regards structure-making/breaking effects, the tra- 
ditional classification of Na + ion as structure-making 
and Cl~ ion as structure-breakin g 23 i 24 ' 44 i 45 appears con- 
firmed in this work. Sodium ions in fact tend to break 
water hydrogen bonds favouring a rearrangement of the 
equilibrium positions of the molecules around the ion 
that results in the formation of definite hydration shells 
(especially the first one). On the other hand, the chlo- 
ride having a greater ionic radius and thus a minor charge 
density, is less effective in breaking the hydrogen bonds. 
Therefore its coordination shells are less defined and more 
disordered as also indicated by the presence of secondary 
peaks in the shell structures. 

We have seen that all peaks intensities of our hydration 
RDFs follow a similar pattern as a function of the density. 
The presence of minima in the peaks heights is probably 
connected to the approach of the system to the spinodal 
line. 

On going from high density toward the spinodal den- 
sity Na + ions partially loose their structure making abil- 
ity while the structure breaking properties of the Cl~ 
ions seems to be reinforced. 

Several extensions of this study can be proposed. In 
particular the analysis of higher densities and concentra- 
tions would be of the uttermost interest in order to in- 
vestigate how thermodynamic and structural properties 
change varying those variables. Indeed varying either the 
concentration or the chemical species of the solutes pro- 
vides a tool to prepare a system with suitable properties, 
for example for cryopreservation. 
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